TO DO: - add plot with eTIV corrections - add CI’s
Let’s load a couple random-ish gamlss models and the dataframe they’re built on
set.seed(12345)
cn_df <- fread("/Users/megardn/Desktop/BGD_Repos/puberty_braincharts/cubic_scripts/bfp_iter/control_filt_og_data.csv", stringsAsFactors = TRUE, na.strings = "")
GMV.int <- readRDS("/Users/megardn/Desktop/BGD_Repos/puberty_braincharts/cubic_scripts/bfp_iter/gamlss_objs/GMV/best_GMV.m3s0_D_mod.rds")
sGMV.re <- gamlss(formula = sGMV ~ bfp(log_age, powers = c(-2, -2, 3), shift = 0, scale = 1) + sex + fs_version + re(random = ~ 1|study) + sex:bfp(log_age, powers = c(1, 2, 2), shift = 0, scale = 1),
sigma.formula = ~ bfp(log_age, powers = c(-2, -2), shift = 0, scale = 1) + sex + re(random = ~ 1|study),
nu.formula = ~ 1, family = GG, data = na.omit(cn_df),
control = gamlss.control(n.cyc = 200), trace = FALSE)
## GAMLSS-RS iteration 1: Global Deviance = 1218118
## GAMLSS-RS iteration 2: Global Deviance = 1218033
## GAMLSS-RS iteration 3: Global Deviance = 1218022
## GAMLSS-RS iteration 4: Global Deviance = 1218017
## GAMLSS-RS iteration 5: Global Deviance = 1218017
## GAMLSS-RS iteration 6: Global Deviance = 1218017
## GAMLSS-RS iteration 7: Global Deviance = 1218017
## GAMLSS-RS iteration 8: Global Deviance = 1218017
## GAMLSS-RS iteration 9: Global Deviance = 1218017
## GAMLSS-RS iteration 10: Global Deviance = 1218017
## GAMLSS-RS iteration 11: Global Deviance = 1218017
## GAMLSS-RS iteration 12: Global Deviance = 1218017
## GAMLSS-RS iteration 13: Global Deviance = 1218017
Call the script holding the plot functions we want to test:
source("plotting_functions.R")
These functions are called from the plotting functions.
sim.data() - takes the dataframe you
built your GAMLSS model on and creates a 2 new dfs with simulated data
(male vs female participants) across the age-range, assigning fs_version
and study values to whatever is most common in the original df. Expects
input df to have log_age, fs_version, and
study. Also preps x-axis labels and defines which centiles
you’ll be plotting.
Returns object sim, which is a list of objects.
sim <- sim.data(cn_df)
names(sim)
## [1] "ageRange" "dataToPredictM" "dataToPredictF"
## [4] "tickMarks_log" "tickLabels_log" "tickMarks_unscaled"
## [7] "tickLabels_unscaled" "desiredCentiles"
centile_predict() - predicts centiles
based on df simulated by sim.data() using the
predictAll() and qGG() functions. Calculates
centiles, 50th centile peak values, and age at peaks separately on male
and female dfs and returns each, as well as an averaged effect across
sexes.
Takes GAMLSS model obj and objects returned from
sim.data(). Returns object pred, which is a
list of objects.
pred <- centile_predict(sGMV.re, sim$dataToPredictM, sim$dataToPredictF, sim$ageRange, sim$desiredCentiles)
## new prediction
## new prediction
## new prediction
## new prediction
names(pred)
## [1] "fanCentiles" "fanCentiles_M" "fanCentiles_F" "peak"
## [5] "peak_age" "peak_M" "peak_age_M" "peak_F"
## [9] "peak_age_F" "M_mu" "M_sigma" "M_nu"
## [13] "F_mu" "F_sigma" "F_nu"
centile_predict.corrected() - same as
centile_predict but corrects out the estimated effects of
fs_version (from mu term) and study (from mu & sigma terms)
parameters.
Expects fs_version to be a fixed effect and
study to be a random effect fit using re()
function!!!
pred.cor <- centile_predict.corrected(sGMV.re, sim$dataToPredictM, sim$dataToPredictF, sim$ageRange, sim$desiredCentiles)
## new prediction
## new prediction
## new prediction
## new prediction
names(pred.cor)
## [1] "fanCentiles" "fanCentiles_M" "fanCentiles_F" "peak"
## [5] "peak_age" "peak_M" "peak_age_M" "peak_F"
## [9] "peak_age_F" "M_mu" "M_sigma" "M_nu"
## [13] "F_mu" "F_sigma" "F_nu"
GGalt.variance - copied from Simon’s Nature paper repo
var <- GGalt.variance(pred$M_mu, pred$M_sigma, pred$M_nu)
makeCentileFan() - basic centile fan
plotting that averages across sex and predicts on Mode(fs_version) and
Mode(study) of original data the gamlss was modeled on. Expects GAMLSS
model, phenotype being modeled, and the name of the original df.
age_transformed parameter set to TRUE or FALSE
makeCentileFan(GMV.int, "GMV", cn_df, TRUE, "sex")
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
makeCentileFan(GMV.int, "GMV", cn_df, FALSE, "sex")
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
makeCentileFan(sGMV.re, "sGMV", cn_df, TRUE, "sex")
## new prediction
## new prediction
## new prediction
## new prediction
makeCentileFan(sGMV.re, "sGMV", cn_df, FALSE, "sex")
## new prediction
## new prediction
## new prediction
## new prediction
makeCentileFan_sex_overlay() - same as
makeCentileFan but with separate centile lines for males
and females
makeCentileFan_sex_overlay(GMV.int, "GMV", cn_df, FALSE, "sex")
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
makeCentileFan_sex_overlay(GMV.int, "GMV", cn_df, TRUE, "sex")
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
makeCentileFan_sex_overlay(sGMV.re, "sGMV", cn_df, FALSE, "sex")
## new prediction
## new prediction
## new prediction
## new prediction
makeCentileFan_sex_overlay(sGMV.re, "sGMV", cn_df, TRUE, "sex")
## new prediction
## new prediction
## new prediction
## new prediction
makeCentileFan.corrected() - centile
fan plot that averages across sexes, correcting for fs_version and study
effects in both centiles and data points plotted.
Expects fs_version to be a fixed effect and
study to be a random effect fit using re()
function!!!
age_transformed parameter set to TRUE or FALSE
makeCentileFan.corrected(sGMV.re, "sGMV", cn_df, TRUE, "sex")
## new prediction
## new prediction
## new prediction
## new prediction
#makeCentileFan(sGMV.re, "sGMV", cn_df, TRUE, "sex")
makeCentileFan.corrected(sGMV.re, "sGMV", cn_df, FALSE, "sex")
## new prediction
## new prediction
## new prediction
## new prediction
makeCentileFan_sex_overlay.corrected()
- same as makeCentileFan.corrected but with separate
centile lines for males and females
makeCentileFan_sex_overlay.corrected(sGMV.re, "sGMV", cn_df, TRUE, "sex")
## new prediction
## new prediction
## new prediction
## new prediction
makeCentileFan_sex_overlay.corrected(sGMV.re, "sGMV", cn_df, FALSE, "sex")
## new prediction
## new prediction
## new prediction
## new prediction
plot.gamlss.var() - plots variance from
predicted GAMLSS model separately for males and females.
age_transformed parameter set to TRUE or FALSE
plot.gamlss.var(sGMV.re, "sGMV", cn_df, TRUE)
## new prediction
## new prediction
## new prediction
## new prediction
plot.gamlss.var(GMV.int, "GMV", cn_df, FALSE)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)